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Abstract 

Weighted essentially non-oscillatory (WENO) and finite volume (FV) methods employ different philosophies 
in their way to perform limiting. We show that a generalized view on limiter functions, which considers a 
two-dimensional, rather than a one-dimensional dependence on the slopes in neighboring cells, allows to write 
WEN03 and 3 rd -order FV schemes in the same fashion. Within this framework, it becomes apparent that the 
classical approach of FV limiters to only consider ratios of the slopes in neighboring cells, is overly restrictive. 
The hope of this new perspective is to establish new connections between WEN03 and FV limiter functions, 
which may give rise to improvements for the limiting behavior in both approaches. 


1 Introduction 

In the finite volume (FV) approach for the numerical solution of hyperbolic conservation laws, flux evaluations 
at the cell interfaces are needed. Since the exact values of the fluxes are unknown at the interfaces, the true flux 
function is approximated by a numerical flux function. Following Godunov [10], this numerical flux function is 
typically evaluated by solving a local Riemann problem at each interface, taking as input the left- and right-sided 
limit of a reconstruction function. If this reconstruction returns simply the adjacent cell mean values, this leads to 
a first order scheme. In order to obtain higher order schemes one can define higher-order reconstruction functions. 
For instance, one can define various types of reconstructions based on only three input values, namely the mean 
values of the cell of interest and its direct neighbors. Evaluating these polynomials at the cell interfaces yields input 
values for the numerical flux function which result in a higher-order accurate scheme. The best order of accuracy 
which can be obtained with only three input values is a quadratic polynomial, resulting in a 3 rd -order-accurate 
linear scheme. Godunov’s Theorem [10] implies that on fixed grids, the high-order approximation of hyperbolic 
conservation laws with linear schemes is impossible without creating spurious oscillations. One way to work 
around this is to limit the order of reconstruction at discontinuous parts of the solution. A large variety of these 
non-linear limiters has been developed to achieve high-order accuracy without creating oscillations [3,11,12,18]. 
The first limiter functions were van Leer’s limiter [27] and Roe’s superbee limiter [22]. These limiter functions 
yield second order accurate reconstructions in smooth parts of the solution and lie in the second order total variation 
diminishing (TVD) region of Sweby [16]. However, due to the TVD property they reduce to l st -order at smooth 
extrema [16]. These classical approaches based on three-cell reconstructions, result in 2 nd -order accuracy [17,28]. 
For smooth solutions however, it is in fact possible to obtain update rules which yield 3 rd -order accuracy. As 
described above, 3 rd -order accuracy is the optimum that can be achieved with only using immediate neighborhood 
information and if only cell averages are stored. Higher orders of accuracy can be achieved only via wider stencils, 
or by storing additional information per cell, as employed for instance in discontinuous Galerkin methods [6,21] 
or jet schemes [5], 

Woodward and Colella [8] proposed a 3 rd -order reconstruction based on a four-point centered stencil. They com¬ 
pute the interface value, and then limit the reconstruction at shocks. Suresh and Huynh [26] go beyond this and 
use a five-point stencil for higher-order reconstruction. The interface values are obtained by limiting a higher- 
order polynomial reconstruction, however, their limiting procedure is costly to apply. Another possibility is the 
use of non-polynomial reconstructions, cf. Marquina [19], who introduced hyperbolic reconstruction schemes, 
and Harten, Engquist, Osher, and Chakravarthy [11], who present essentially non-oscillatory (ENO) schemes. The 
idea behind ENO schemes is to divide the stencil of interest in smaller sub-stencils with piecewise polynomial re¬ 
constructions on each sub-stencil. The scheme then applies an optimal stencil selection procedure, which chooses 
the locally smoothest stencil for the reconstruction. The scheme avoids interpolation across discontinuities, how¬ 
ever, the final approximation does not contain all available data. To overcome this drawback, weighted essentially 
non-oscillatory (WENO) schemes were developed by Liu, Osher, and Chan [18], Weighted ENO schemes use a 
convex combination of all candidate stencils to obtain the reconstruction of the interface values. Later, Jiang and 
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Shu [12] further modified and improved the WENO scheme by proposing a new way of measuring the smoothness 
of a numerical solution and thereby increasing the order of accuracy. 

Artebrant and Schroll [3] developed a local double-logarithmic reconstruction which only needs three input values. 
Based on their work, Cada and Torrilhon [4] introduced a limiter function, which yields 3 rd order accuracy in 
smooth parts of the solution and at extrema. The authors also introduce a criterion to distinguish between smooth 
extrema and discontinuities in order to avoid extrema clipping. This work has been further investigated in [23],in 
order to find a parameter-free smoothness criterion. 

This paper continues their work on determining a parameter-free smoothness criterion and incorporating it in the 
reconstruction procedure. We present a new 3 rd -order limiter function in the FV setting and then compare it to 
different 3 rd -order WENO schemes with the goal that this contributes to a better understanding of the nature of 
non-linear schemes. 

The time discretization of all schemes is implemented using the 3 rd -order TVD Runge-Kutta method developed by 
Shu and Osher [24], 

The paper is structured as follows. In Section 2 we review the basic formulation of the problems in question 
and interface reconstruction in general. In Section 3, we discuss FV reconstruction methods containing limiter 
functions and we provide a new smoothness indicator. Section 4 recalls the WENO scheme as first introduced by 
Liu et al. [18] and later improved by Jiang and Shu [12], We focus on WENO schemes which only need three 
points for each reconstruction, as from now called WEN03. Section 5 places the introduced methods into a uni¬ 
fying setting and compares the novel limiter function with different WENO schemes. Finally, Section 6 presents 
some numerical results which demonstrate the potential of the proposed schemes and conclusions are given in 
Section 7. 

2 Basic Formulation 

We are interested in the numerical approximation of a Cauchy Problem of the hyperbolic conservation law 

d r u(x,t ) +d x f(u(x,t)) = 0, (2.1a) 

m(x,0) = uq(x), xGK. (2.1b) 

in one space dimension, equipped with the initial condition uq(x). where u = (mi, , u s ) T and the Jacobian matrix 
A (u) = has ,v real eigenvalues. For the sake of simplicity, we restrict our discussion and analysis to the scalar 
case 5=1. However, the developed ideas are also applicable to systems (s > 1), in the same way other limiters 
extend from s = 1 to s > 1. We consider a regular grid in space, with the positions of the cell centers denoted by 
Xi, i € Z and with uniform space intervals of size Ax. The grid cells are defined by C, = [*/-i/2> V+ 1 / 2 L where 
Xj±j = Xi ± jAx. Finite volume (FV) methods aim at approximating the cell averages of the true solution of (2.1) 
with high accuracy, see e.g. [17]. The cell average of the true solution «(•,•) is given by 

J rx+Ax/2 

U(x,t) = — / u(s,t)ds. (2.2) 

Ax Jx—Ax/2 

With this definition we call U" = U(x/.t") the cell average of the true solution u in cell C, at time t n . The goal is to 
find an update rule to advance approximate cell averages from a given time t n to a new time t" : 1 = t“ + At, such 
that the true cell averages are approximated with high order of accuracy. In addition, the approximate solution 
should not develop any (relevant) spurious oscillations. Integrating Eq. (2.1a) over the cell C, and dividing by Ax- 
yields 

^ (2.3) 

which is still exact. We now want to find a solution approximation m" satisfying i/" « U". The quality of the 
approximation m, depends on the accurate approximation of the fluxes at the cell boundaries f(ii(x [+ \ n-t)). This 
is achieved by constructing a numerical flux function fiu. : v) that is Lipschitz continuous and consistent with the 
(true) flux function, i.e. /(m,m) = f(u). The numerical fluxes at the boundaries of cell C, are then given by 

fi+ 1/2 = /(“i+l/2’“/+l/ 2 )> 

/U/2=M ( :| /2 ,«H /2 )’ 
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(2.4a) 

(2.4b) 
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(a) From the three cell averages (black horizontal 
lines), functions are reconstructed (red lines, a 
linear function in this case), these functions are 
then evaluated at the left- and right-sided cell 
boundaries. 
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(b) The normalized gradients are computed as the 
differences between neighboring cell averages, 

8-, i = u, + 1 — Uj and 8- i = «; — iij_ i, 

' ' 2 1 2 
respectively. 


Figure 1 - Basic setting for the reconstruction of the interface values on a 3-point-stencil. 


where u ; 2 and | , 2 are approximations to the solution values at the cell boundary x J±1/ / 2 , left-sided (( )) and 
right-sided (/ 1 '), respectively, see Fig. 1. The evolution of cell averages is thus given by 


dui 

dt 


~Av (^Tl/2 fi—l/l) ■ 


(2.5) 


The accurate reconstruction of these left and right interface values at the cell boundaries *,± 1 / 2 , see Fig. 1, is the 
crucial point in this process. We are particularly interested in numerical schemes that find the approximate solution 
values u\ ' j, 2 and , 2 that correspond to cell C, by using only information of the cell C, and its immediate 
neighbor cells C,_i and Q+i. The restriction to immediate neighbor cells provides local update rules. This locality 
is, among other advantages, beneficial for low-communication parallelization, and ensures that few ghost cells 
must be provided near boundaries. 

The key ingredient to getting 3 rd -order accuracy is the way of reconstructing function values at cell boundaries 
x j± 1 based on cell averages. The reconstructed values u j± 1 are then provided as input values for the numerical flux 

function /(•,•) in Eq. (2.5), following the standard FV methodology [16]. The focus of this paper is solely on the 
actual reconstruction of the solution u j± 1 at the cell interfaces. For the sake of simplicity, we shall drop the 
Since we consider three cells for the reconstruction, we define the left and right interface values to be determined 
by functions L and R: 


u [ I =L(u i - l ,u i ,u i+ i) 
l ' 2 

u [+ \ =R(ui-i,Ui,u i+ i). 


(2.6a) 

(2.6b) 


3 Finite Volume 3 rd order Limiting 


In this section, we review limiting in FV methods using the three-point stencil for the reconstruction of and 

ilj ± j , 2 . For the computation of the numerical flux functions, the approximate solution at the cell interfaces x ( - ± | / 2 
is needed. For the cell x,. we define the left and right interface values by 


“z+i /2 = «i+ \ 0(0/)5,- + 1 , (3.1a) 

“l-i/2 = «/ - \ (3- lb ) 
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respectively. Here, 0 is a univariate, non-linear limiter function depending on the local smoothness measure 
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(3.2a) 
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where 5 ;+ r = w, + i — m,. 

(3.2b) 

and S. 1 = tij — i7,-_ 1 
' 2 

(3.2c) 


are the differences between neighboring cell averages, cf. Fig. 1. Limiter functions switch the reconstruction 
to high-order accuracy in smooth parts of the solution and to lower-order reconstructions near discontinuities 
[16], This means, limiter functions contain some sort of smoothness measure. In Eq. (3.1), the choice of 0(0,-) 
determines the order of accuracy of the reconstruction and therefore of the resulting scheme. 

Property 3.1. General Properties of Limiter Functions 

(i) If(j) passes continuously through 0 = 1 with 0 ( 1 ) = 1 , then the resulting scheme is at least 2 nd order accurate 
in sufficiently smooth, monotonous regions of the solution. 

(ii) If 0 satisfies the conditions 

0 < 0(0) < max(O,min(2,20)) , (3.3) 

the numerical scheme is total variation diminishing (TVD), and thus does not create spurious oscillations. 


Proof. See e.g. [17]. 


□ 


There is a variety of schemes which reconstruct based on three points and obtain 2 nd -order accuracy. These are 
the classical schemes in the Sweby setting, such as Superbee, van Leer and others [16], These schemes use the 
information of the three cells to compute a linear reconstruction function, see e.g. [28], Indeed, the full 2 nd - 
order reconstruction u\ + \p = if, + ( H ' + 1 oa ^'~ 1 ) can be rewritten in form of Eq. (3.1) with the limiter function 

0(0) = 4+ One can easily check that this limiter function satisfies the property 0(0 1 ): 0 1 0 (0) and therefore, 
for this limiter, Eq. (3.1) can be rewritten in the formulation 


(-) - Ax 

U i+ 1/2 - + y 


..(+) 


Ax 


* 1 - 1/2 — 2 ° 


(3.4a) 

(3.4b) 


with the right-sided slope <7, = 0(0,)5 + i/A y, see e.g. [17]. Even though this formulation is widely used, we will 
continue with the formulations (2.6) and (3.1), as introduced in [4,9]. The aim of this work is to discuss schemes 
which use the three-point stencil to achieve 3 rd order accurate reconstructions of the cell-interface values. One 
possibility is to construct a quadratic polynomial pfx) in each cell Q. Applying the computed polynomial to 
x,-±i/ 2 , we obtain the interface values 


ll \+\/2 = P’( x i+ 1 / 2 ) = 3 M '+i + ^ gKi-i (3-5a) 

M !-l /2 = Pi{ x i- 1 / 2 ) = \ U '~ l + \ Ui ~ j+‘+l- (3.5b) 

It turns out that these expressions can be written in the form (3.1) to obtain the non-limited 3 rd -order reconstruction 


03 (ft) := 


2 + 0 / 
3 


(3.6) 


with 0, given by Eq. (3.2a). This formulation results in a full-3 rd -order-accurate scheme for smooth solutions, 
however, causes oscillations near discontinuities. This can be seen either by noticing that 03(0,) does not lie in 
the TVD region, see Fig. 2, or as a direct consequence of Godunov’s Theorem [10], since a linear scheme of more 
than 1 st order cannot be monotone. Since oscillations should be avoided, limiter functions have been introduced, 
which apply the full 3 rd -order reconstruction (3.6) at smooth parts of the solution and switch to a lower-order 
reconstruction close to large gradients, shocks, and discontinuities. In [3], Artebrant and Schroll present a limiter 
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Figure 2 - Logarithmic limiter function 0 as with < 7 = 1.4 (black solid line), its approximation 0c t (blue 

solid line), the TVD-version of 0 ct which does not have the non-zero part for 0 < 0, denoted by 
0ct, tvd ( re d solid line with red circles) and the full-3 ld -order reconstruction 03 (black dashed 
line). 



function, which can be formulated as 0as (0<,<?), see [4], based on a local-double-logarithmic reconstruction. This 
function does not solely depend on 0 , but also contains a parameter <7 which significantly changes the reconstruction 
function. The authors recommend < 7 = 1.4 and demonstrate that for q —» 0, the logarithmic limiter function reduces 
to 03 (0,). Their limiter function reads 


0AS (9i,q) 


2 p[(p 2 - 2 p 0 / + 1 ) log(p) - (1 

(p 2 -l)(p- 1) 2 


p=p(9i,q) 


2 

1 + | 0 *| 2 «' 


mp 2 


i)] 


The downside of 0 as( 0!,<?) is its complexity, which renders the evaluation in each cell expensive. Cada and Tor- 
rilhon [4] derive, in an ad-hoc fashion, a limiter function 0ct(0<) which is based on 0 as- This function overcomes 
the drawbacks by approximating the properties of 0 as and reducing the computational cost. It reads 


0ct( 0;) = max 


^ 0 ,min 


^ 0 3 ( 0 ,),max 



0 <,min (2 01 , 03 ( 0 ,), 1 . 6 ) 


(3.7) 


and is shown in Fig. 2 together with 0as(0<, 1.4) and 03 ( 0 ,). Note that (pa does not lie within the strict TVD 
bounds, i.e. it breaks with 3.1 (ii). The TVD property 3.1 (ii) can be achieved by considering only those parts of 
0CT where 0 > 0, and setting it to 0 elsewhere, leading to the limiter function 


0CT,tvd(0) =max(O,min(20,0 3 (0),1.6)) . (3.8) 

However, the motivation for keeping the non-zero part in the construction of 0ct(0) f° r 0 £ [—2,0] is to avoid the 
so-called extrema clipping. This is the effect occurring close to minima and maxima, where the normalized slopes 
5,-± 1/2 are of the same order of magnitude but have opposite signs, i.e. 0,- « —1. In this case, classical limiter 
functions that fully lie in the strict TVD bounds yield zero and thus generate a l st -order accurate scheme. This 
undesirable reduction in accuracy is avoided when the non-zero part of 0 is included. Therefore, (per possesses 
better smoothness properties near 0 = — 1 than 0ct, tvd since here, (pa (0 ) = 0 3 (0) but 0ct, tvd ( 0) = 0. For more 
details, see [4], where the non-zero part was first introduced. 

Nonetheless, it might occur that the discretization of an extremum is such that one of the consecutive slopes is 
approximately zero. In this case, 0 —> 0 or 0 —► ±°° and the interface values wy- 1/2 are again approximated by 
the cell mean values which yields a l st -order scheme. That is, a zero-slope is interpreted as the onset of a 
discontinuity, even though it might in fact be the magnified view of a smooth extremum. This undesired case 
demonstrates that a criterion is needed which can differentiate a smooth extremum from a discontinuity or steep 
gradient. In the framework considered in this paper, the criterion should only depend on the available information 
of the compact three-point stencil. Furthermore, it has to detect cases when switching to the 3 rd order reconstruction 
is safe, even though one of the normalized slopes is zero. We assume that using the 3 rd order reconstruction is safe 
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(a) Limiter function Hex tvd (<5,_ i ■ 8 i+ 1) satisfying 
the strict TVD bounds. 


(b) Extended limiter function Hqt (<5-_ i. 8 i+ 1 ). 


Figure 3 - Different limiter functions in the two-parameter framework. 


if the non-zero slope is ’small’. In turn, if the non-zero slope is not small, we assume to be near a discontinuity or 
large gradient, and the order should be reduced. The main focus of Section 3.1 is to determine what ’small’ means 
and to define a suitable smoothness indicator r|. 


3.1 Interpretation in 2D Slope Domain 

From the discussion above, it is clear that such a switch function r| has to explicitly depend on both normalized 
slopes 8j z 1 - 1 / 2 , Eq. (3.2b, 3.2c). The classical approach of only considering the ratio 0,-, Eq. (3.2a), of neighboring 
slopes is overly restrictive because part of the information (the actual magnitude of the two slopes) is discarded. 
This is why we reformulate all limiter functions 0 in a two-parameter framework and obtain the new formulation 
for the reconstructed interface values (see Eq. (3.1)) as 

“U/2 = ««+ 5 H ($_i, 5,+1). (3-9a) 

»£\,2=*i-\m + \A-\) (3.9b) 


with the limiter function H explicitly depending on both normalized slopes. The old limiter function (p ( 0,) can of 
course be rewritten in the new form of the two-parameter function H by setting 


H(8 iA ,8 i+ i):=m)8 i+ i 



(3.10) 


In this setting, the full-3 rd -order reconstruction 03(0,) = (2 + 0,)/3, given by Eq. (3.6), now reads 


H 3 (8._ 


A + 0 = 


25 ;+ i + 8 ; 


(3.11) 


This formulation has the advantage that there is no division by the normalized slope 5 I± | ; / 2 . Thus, a possible 
division by a number close to zero is avoided. Fig. 3a shows the limiter function //ct.tvd which satisfies the strict 
TVD bounds. This can clearly be seen by the zero parts for sgn( <5 i) ^ sgn( 8 j+ 1 ). Fig. 3b shows the extended 
version I Icy in the two-parameter setting. On the coordinate axis where 8._ i = 0, i.e. 0,- = 0, the limiter function 
I Icy returns zero, meaning that it yields a l st -order method. The same holds for the coordinate axis where 5 (+ 1 = 0. 
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For two consecutive slopes of approximately the same order of magnitude, i.e. around the diagonals, Hqt recovers 
the 3 rd -order reconstruction . This is the case for 8 _ 1 « 8 ,1 as well as for 8 1 ss 5, 1 , contrary to flcTTVD 

which returns //; only in the latter case and 0 for 8 1 sa — 8 ., 1 . 

I_ 2 ! +2 

The limiter function presented in [4] is not symmetric with respect to the diagonals. This means that in some cases 
Hcj{8\ , Si) = /Tt((5|, 82 ) but Hct{—82, — 81 ) f H^—Si, <5|), cf. Fig. 4. This asymmetry is not natural because 
these two situations are equally smooth or non-smooth. We therefore correct this feature by defining the following 
3 rd -order limiter function fluff). 1 , <Y ( 1 ) 

// 3 l (<5._i,c5. + i) =sgn(<5 ;+ i) max(0 , min(sgn(5. + 1 ) #3 , max(— sgn(5 ;+ 1)1 , (3.12) 

min (2 sgn (5 j+ 1 ) 5._ 1 , sgn (5. + 1 ) ,1.515,. + 1 1)))) 

This new limiter function treats symmetric situations in the same manner, i.e. if HmJ 8 \ , 82 ) = /V 3 (<5 1 , 82 ) then also 
//ii f-cT. — 5]) = H$(— 82 , — 5i), cf. Fig. 4. The difference between the functions //f i and W 31 can also be seen in 
Fig. 7. “ 


Expressing the interface values in the more general form (2.6), we can determine some properties for L(-, ■, •) and 
/?(-, -, •), as done in [9]. These properties are valid for all limiter functions presented so far, i.e. for //;. Her, //;m.. 


Property 3.2. Homogeneity 

Multiplying the arguments ofL and R in Eq. (2.6) by the same real number A multiplies the interface values u. | 


<+1 


and u h j, respectively, by the same constant A. This is, L and R are called homogeneous, i.e. linear along each 


i~\~ 


line through the origin in the (5._i , 8 . 1 ) plane. 

I -1 l~\~ ri ' 


7(Am,Av, Aw) = A J(u,v,w), J £ {L, R }, A £ K. 


(3.13) 


Property 3.3. Translational invariance 

Adding a constant A to the arguments of (2.6) adds the same constant A to the interface values. This is, L and R 
are called translationally invariant. 

/(m 4 -A,v + A,w + A) =/(m, v,w) + A, Jg{L,R}, A £ K (3.14) 


Property 3.4. Left-Right symmetry 

Exchanging the first and third argument of (2.6) interchanges the left and right interface values, (cf. [9] for more 
details and figures) 


R(w,v,u) = L(u,v,w). 


(3.15) 


u 



(a) This situation is treated as a possible discontinuity: 
Hct f Hi ■ 


U 



(b) This situation is classified as smooth: Hqt = H 3 


Figure 4 - Two situations which are reflections of each other, treated differently by flier. 
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-g- smooth solution 
cell averages 

-a- discontinuous solution 


Figure 5 - Identical cell mean values, corresponding to a smooth solution, or to a discontinuous solution. 


. such that 


Lemma 3.1. 

If properties 3.2 to 3.4 are satisfied, there exists an appropriate limiter function if/ : 

. 1 

L(u,v,w)=v+ -iff 

, , 1 
R(u,v,w) =v- -iff 

Lemma 3.2. With properties 3.2 to 3.4, it is easy to verify that 

(i) computing the interface values j /2 with L and R, given by Eq. (2.6), or in the form of Eq. (3.1), are 
equivalent. 

(ii) For S i+ i 0 the formulations of the cell interfaces in the one-parameter framework, Eq. (3.1) and in the 


u — V > 

) (w - v) 

(3.16) 

W — Vj 

w — V ^ 

) (v-k). 

(3.17) 

U — Vf 


two-parameter framework, Eq. (3.9), are equivalent. 

Proof. We will only show (i) for w| + jthe other cases are similar, 
(i) Setting ( u,v,w ) = yields 


1 


1 


(2-6) , , _ _ _ n (3.1) _ , V Ui-i-Ui\ _ _ 1 

= L{ui-i,Ui,u i+ i) = Uj + -\f/\ —— ^-\{ui+i-Ui)=Ui+-0(8i)b j+ i 


V. . , ,r. 

+ / ' ' '■£ ' \Ui+l-Uiy 

where <5 ;+ 1 = ti, + 1 — Cij and y/(iii-i —Uj/uj + \ —Hi) = y/(— 0,) =: </>(0i) applying Eq. (3.2a). 
(ii) Due to the homogeneity of L and R, we can easily see with Eq. (3.9) that 

5, i 

H(8 i _i,8.,i) =H(—^,1) 8- i =H(G h 1) 5 = m) 5 i+ i- 




(3.18) 


□: 


3.2 A New Smoothness Indicator 


Given three cell averages, it is in general not possible to determine whether these points represent the onset of 
a discontinuity or the magnified view of an extremum. An example is shown in Fig. 5 where one set of cell 
mean values could be obtained by the smooth function (green, square markers) or the discontinuous function (red, 
diamond markers). As stated in Sec. 3.1, the two-parameter setting is the necessary prerequisite for the definition 
of such a criterion. Cada and Torrilhon [4] proposed the function 


tJct(5 ; _i ,5 i+ i) 



(rAr) 2 


(3.19) 
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Figure 6 - Double logarithmic plot of the L\ - and L^-errors versus the number of grid cells of the solution 

(c) 

obtained with © r| advected until f en( j = 1 with CFL number v = 0.9. The numerical solutions 
are shown for different values of r, see Eq. (3.19). Left: L\ -error, right: -error. 


This switch function defines an asymptotic region of radius r around the origin in the (5. >, 5 i )-plane. Within 
this region the limiter switches to the 3 rd -order reconstruction. 

The authors of [4] then modified the structure of the limiter function 0 ct to include the asymptotic region around 
the origin. The combination, denoted by the superscript (c), is defined as 



if r/cT < 1 
if T ) ct > 1 ■ 


(3.20) 


It is possible to make this function Lipschitz continuous, by introducing a small transition region and a linear 
function, cf. [4] for more details. This limiter function, combining (pa and the switch function r/ci, has been 
successfully employed in e.g. [13,14,20]. In [4], the authors did not provide a general formulation of the param¬ 
eter r which determines the size of the asymptotic region. Instead it was chosen ad hoc, in a problem-specific 
way. To obtain some generic idea about suitable choices for r, we conduct a numerical test, applying (j) 1 ^ to the 
advection equation u t +u x = 0 with smooth initial condition uq(x) = sin( 7 tv) for different values of r. Fig. 6 shows 
the double logarithmic plot of the L\ - and -errors versus the number of grid cells of the solution advected until 
tend = 1 with CFL number v = 0.9. For this smooth test case, we see that larger values of r, corresponding to larger 
asymptotic regions are favorable. For smooth solution this makes sense because increasing values of r corresponds 
to increasing the region of directly applying the full-3 rd order reconstruction (p 3 . From Fig. 6 we can deduce that 
for smaller r, a finer space discretization is needed to obtain 3 rd order accuracy. 

Motivated by the limiter function of [4], we want to define a new smoothness indicator without artificial parame¬ 
ters. As already sketched in [23], a promising potential to distinguish discontinuities from smooth extrema is by 
measuring the magnitude of the vector (5 ( . _i,8 i+ i ). If this vector is sufficiently small in some appropriate norm, 

the reconstruction is switched to the full-3 rd order reconstruction, even if one of the lateral derivatives may be 
vanishing. 

Lemma 3.3. In the vicinity of an extremum go, for \x i-So |< Ax', the following relations hold for each time t": 



2 

1 


< s/c max \u"{x,t n )\Ax 2 with c = - + &(Ax), 

xen\Q. d 2 

< c max | u" (x, t n ) | Ax 2 with c = 2 + 6(Ax). 


(3.21a) 

(3.21b) 


Here, fl is the computational domain, and is a set of points where the solution is discontinuous. 
Proof Let us recall the definition of U(xf), Eq. (2.2), 

J rx+Ax/2 

U(xf) = — / u(s,t)ds , 

Ax Jx—Ax/2 


9 












where u is the exact solution. The following properties hold for U(x,t) 


U 1 (x,t) 
U"(x,t) 


u(x + ^,t) — m(x — 4^,f) 
Ax 

u'(x+ ^f,t) — u'(x— ^r,t) 
Ax 

up = u( Xi ,t n ), 


(3.22a) 

(3.22b) 

(3.22c) 


where U ( Xj,t n ) is the cell average of the exact solution at time t" in cell C, with cell center x,. 

Eq. (3.21a) can be proven regarding the following formulation of definitions (3.2b), (3.2c), with a constant a which 
is going to be specified. For the sake of simplicity we shall neglect t n for the rest of the proof. 


( U (xi + Ax) — U (x /)) 2 + (U (xj) — U (x,- — Ax )) 2 
(a Ax 2 ) 2 

1 / U (xj + Ax) - 2U(xi) + U (xj - Ax) \ 2 2 / f/(x,- + Ax) — U(xj) \ /f7(x,) — U(xj — Ax) 

a 2 \ Ax 2 ) a 2 Ax 2 \ Ax ) \ Ax 


(3.23a) 

(3.23b) 


a Taylor expansion of the functions U (x,- ± Ax) around x, yields 


1 ( U"{xj ) \ 2 2 ( U'(xi ) \ 2 2 U'{xj)U'"(xj) 

2 \ a ) Ax 2 \ a ) 3 a 2 


6 (Ax 2 ). 


In the vicinity of an extremum |o, for |x,- — Co < Ax, the derivative satisfies 

U'{ Xi ) = £/'(§ 0 ) + £/"(&)(*,- - §,) + 0{{xi - ^ 0 ) 2 ) 

= U"(^o)(xi — |o) + ^(Ax 2 ) 
^|t7'(x ! -)|<|t7"(^o)|Ax+^(Ax 2 ). 


Since we are interested in ({/'(x,)) 2 , we find that 


(t/'(x,)) 2 < (t7"(^)) 2 Ax 2 + ^(Ax 3 ). 


(3.23c) 


(3.24a) 

(3.24b) 

(3.24c) 


(3.24d) 


Therefore, Eq. (3.23) reduces to 


(t/(x,- + Ax) — f/(x;)) 2 + {U(xj) — U(Xi — Ax)) 2 1 fu"{xi) 

(aAx 2 ) 2 “2^ a 



2 

+ &(Ax). 


(3.24e) 


Consider the computational domain £1, and the set of points £i,/, where the solution is discontinuous. Setting 
a = maXj.g^o, \U"(xi,t”)\ with the cell centers x,-, leads to 


((7(x,- + Ay) - U(xj) ) 2 + {U (x,-) - 6' (xj - Ay) ) 2 < 5 + 
(max^. \U"(x,t n )\Ax 2 ) 2 2 


(3.25) 


Since the numerical solution u" is a 3 rd -order-accurate approximation of the true solution, i.e. u" = U(xj.t n ) + 
i^(Ax 3 ) and 


( 3 . 22 ) 

a =max\U"(xi,t n )\ < max max \u"(^,t)\ 

Xi 

holds true, this shows Eq. (3.21a). In a similar manner, Eq. (3.21b) can be proven. □ 

Remark 3.1. Often, the exact value o/max* \u"{x,t n )\ is not known or it is too expensive to compute. In these 
cases, a different estimator needs to be found. 
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(a) In many applications, one has some estimate of the largest second derivative of the solution, even if one does 
not know the solution itself. 

(b) If a good estimate of the solution at time t n is unavailable, one can use the initial conditions uq (x) as 
an approximation of a. Note that for a conservation law of the form u t + f(u) x = 0, certain information 
about u xx can in fact be inferred from the initial data (where the solution is smooth). For instance, second 
derivatives are actually preserved at extremal points. To see this, consider the equations that u x and u xx 
satisfy, namely (u x ) t + f (u)(u x ) x = - f"(u)(u x ) 2 and (u xx ) t +f(u)(u xx ) x = -f'"(u)u\ - 3 f"(u)u x u xx . In 
both equations, the left hand side represents the convective derivative along characteristics. Therefore, 
extremal points (u x = 0) are just moved with the characteristics, and u xx remains unchanged from its value 
at initial time. 

From now on, we will use a = max ve Q\n rf \ u o Ml• 

Lemma 3.3 makes a statement on the magnitude of the normalized slopes. Note that the upper bound only de¬ 
pends on the grid size Ax and the initial condition uq. From this lemma we can now deduce the definition of the 
smoothness indicator. 

Definition 3.1. The switch function rj which marks the limit between smooth extrema and discontinuities is defined 
by 


with 






(3.26) 


a= max |ma(x)| 
xeQ.\Q. d 


(3.27) 


FI and Fid defined as in Lemma 3.3. 

Note that a is proportional to 1 /Ax 2 and therefore, I] is a non-dimensional quantity, so that the overall scheme is 
not affected by changes of units. 

The idea of using the largest second derivative of the initial conditions to relax limiting near extrema has already 
been proposed by Cockburn and Shu [7] in the context of discontinuous Galerkin methods. They use the constant 
M 2 = max v if/(x) to overcome the degeneration to first order at critical points: Mi is used to estimate the mag¬ 
nitude of the reconstructions / 2 i M ;±i/ 2 - ^ the reconstruction is smaller than a certain bound, the high order 
reconstruction is used, otherwise, a limiter function is applied. Note that in [7] the switch function is based on the 
reconstructed values i/ +| */ 2 > M j±j/ 2 ’ whereas in this work, we consider the normalized slopes 5 i±1 / 2 . 


With Def. 3.1, Lemma 3.3 states that in the vicinity of smooth extrema, ij < 1 holds. Combining this information 
with the new limiter function // 3 L, we use this result to define the combined limiter, denoted by the superscript (c), 


<A- 


A + i):= 


| w 3l(5. 


■Vi) 


if 17 < 1 
if 77 > 1 . 


(3.28) 


The resulting two-dimensional limiter function for this approach is shown in Fig. 7b. Note that in the same manner 
as for Hf//, Lipschitz-continuity of the combined limiter II < // [ l (8 j 1 , 8 j+ 1 ), Eq. (3.28), can be achieved by using a 
continuous switch function for the transition between the limited and non-limited reconstruction. 

Remark 3.2. The proposed limiter H^(8 i _i,8 j+ i) still satisfies properties 3.2 to 3.4 , i.e. it is homogeneous, 
translationally invariant and fulfills the left-right symmetry. For the homogeneity, this can be seen by the fact that 
the stretched vector (kui \,ku,,kui + \) leads to the undivided differences (k<5_i ,k 8 j+ 1 ), which are stretched by 
the same factor. Also, obsen’e how (kui-\,kui,kui+\) affects a, Eq. (3.27), 


a(kuo) = max \ku'/(x)\ = |k| max |mq(x)| = |&|a(no)- 

xeQ\Q d xeQ\Q d 


(3.29a) 
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(a) Combined limiter function Hq£, [4], (b) New combined limiter function . 

Figure 7 - Limiter functions combined with the full-3 rd order reconstruction in a region around the origin. 


With definition (3.26), we obtain that the decision criterion is homogeneous, i.e., 

V(k8 iA ,k8 i+l2 ) = r 1 (8 iA ,8 i+li ). (3.29b) 

Using the fact that Ht, (5._ i , 5., i ) and 8 . i) are homogeneous functions, finally leads to 

l J l~\~ r\ l r\ 1 I ^ 


ft &i_ i , k 8 i+ 1 ) — k H^l (S t _ i , 5 j+ 1). 


(c), 


(3.29c) 


For the translational invariance, note that the shifted vector (k + w,_ i , k +ut,k + u l+ \) leads to the original for¬ 
mulation (5- i, 8 j+ 1 ) and does not affect a either. Contrary to H^f, the limiter function _ i, 8 j+ i ) is not 

homogeneous, because ricrikSj i,k 8 j+ 1 ) = k 2 t7cr(£> ; _ i , <5 (+ 1 ) Tjcr(5 i _ i , 8 i+ 1 ), compare to Eq. (3.29b). 


4 Third order WENO 

In this section we briefly review the 3 ld -order weighted essentially non-oscillatory (WENO) scheme for one¬ 
dimensional scalar conservation laws. We consider WENO schemes because they represent a very popular method 
for approximating the solution of hyperbolic conservation laws. Since we want to compare WENO with the limiter 
functions introduced in Section 3, we focus on WEN03. This is the family of schemes which only use three points 
for the reconstruction of the left and right cell-interfaces id j^ 2 and w| + j , 1 . WEN03 reconstructs an estimate of 

the solution u(x. , i) from the cell averages and ti,_i. The procedure, introduced by Jiang and Shu [12], 

'+2 

is to start just as for the r-th-ordcr ENO scheme [11], where r = 2. This is, we consider the two-point stencils 
So = Q-i U Cj and .S'] — (.] I—J Ci j i and define on each stencil 2 -order accurate approximations. 

• po(x i+ i) = —\ui -1 + 2 Uj , where po is the linear interpolation of (jti_ i, m,-_ i ) and and 

• P 1 ( x i + 1 ) = 3 m / + h u i+ 1, where p\ is the linear interpolation of and (x i+ i,M,- + i). 

For the sake of simplicity, we only consider the procedure for the right cell interface x-, i . The left interface follows 

i+i 

in a similar manner. 

Based on an r-th order ENO scheme, the best one can get is a (2r — l)-th order WENO scheme, i.e. in our case 3 rd 
order [12], To obtain this, the WEN03 estimate of u (x. + 1 ), called u j+ 1 , is defined as a convex combination of the 
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(4.1) 


two 2 nd -order estimates po (x._ 1 ) and p\ (x._ 1 ): 

“i+> = ) +w,. + ipi(* j+ i). 

The weights satisfy w ; . ± i > 0 and w._ i + vv . + 1 = 1. There is a particular choice of weights that generates a 3 rd - 
order-accurate approximation to u(x j+ i ), namely vv ( _^ = Y t _i = 3 , and w j+ i = tj. + i = |. This approximation is 
obtained when interpolating (x,-, w,-), and (x i -+i,m i -+i) by a quadratic polynomial, and evaluating it at 

x,., 1 , see Eq. (3.5a). 

The philosophy of WENO is to design formulas for the weights w. ± 1 , such that in smooth regions, one has w j± 1 ~ 
y j± 1 , while in regions of large gradients, more weight is given to the approximation that generates fewer spurious 
oscillations. This is achieved using smoothness indicators pj±in- They are defined using the normalized slopes 
between neighboring cells 5., 1 : 

The final weights w ,, 1 are defined by 


where 


In equation (4.3), there are two parameters which need to be further detailed. The integer p € N, which Liu et 
at. [18] suggest to set p = r, the order of the base ENO scheme. Jiang and Shu [12] claim that for r = 2,3, setting 
p to 2 is adequate. In this work, we will set p = 2. 

The other parameter in Eq. (4.3) is e, a small positive number, originally introduced to avoid the division by 
zero [12]. We suggest that rather than fixing £ to some constant, it should depend on the spacial discretization Ax. 
This will be discussed in more detail in Sec. 4.2. 


= ($**)■ 


«/±i 

W., 1 =- - - 

l± 2 a. i+a,.,i 




a,± 2 (e + P i+i )P 


(4.2) 


(4.3) 


4.1 Interpretation of WEN03 in 2D Slope Domain 

A brief calculation yields for the WEN03 reconstruction 




= w 

l 


) +vv /+ 

\P^ x i+\) 


(■ 

-5“f- 

1 + \ui 

)+ w i+\ (2 


1 

/ 

(Uj — Ui 

-1 )+ w i+i( 

= Ui 

+ 2 

( W i-i 

= Ui 

1 

+ 2 


8 i ,+ 

w i+\ 8 i+\) 

only 

depend on the slopes 5 ;± 1 , t 


uuivv uiv vy rv • I 1 uiviiiov. 

. 2 

form (3.9) with the limiter function W W r,N 03 . whose particular form depends on the parameters p and e. Explicitly 
written, it is 

„ /s s x 

^WEN03 (5-_ 1 ,5. 1 ) = - t ; C- - NT' . /c- --- ■ (4-4) 


(e + (5,._i ) 2 )- p + i(e + (8 i+1 ) 2 )-P 


In the vicinity of (0,0), i.e. for |5, ± 1 1 <C £, one has in leading order that 


H«(8 t i,8 i+ i) - 31 + l§ i+ i . 


(4.5) 


As mentioned before, this linear function is homogeneous and results in a 3 rd -order-accurate approximation. On 
the other side of the spectrum, if |5 ;± 1 | £, one has in leading order that 




(4.6) 


13 





Figure 8 - Double logarithmic plot of the L\ - and L x -errors versus number of grid cells of the solution 

obtained with //weno 3 advected until f en( i = 1 with CFL number v = 0.9. Numerical solution for 
different values of £, see Eq. (4.4) and (4.3). Left: Li-error, right: L^-error. 


This function is also homogeneous, i.e. Hj > (k8 i i./rcCi) = kH s> (8 i _ i, S j+ 1 ). That means, along each line 
through the origin in the (<5-_ i, 8. + i) plane, it is linear, thus resembling the behavior of traditional FV limiters. 

4.2 WENO Smoothness Indicators 

The choice of the weights which depend on a, ±\n, is crucial for the order of accuracy of the resulting 

scheme. Its precise value influences the behavior of the limiting when j3 i± t/2 is close to zero. Clearly, this is 
of particular importance near extremal points of the solution. To point out the importance of £, we repeated the 
numerical experiment conducted in Sec. 3.2 with the 3 rd -order WENO scheme. This is, we apply //\vh \03 to the 
advection equation u, + u x = 0 with smooth initial condition iiu(x) = sin( tta - ) for different values of £. The result 
for fend = 1 and CFL number v = 0.9 is depicted in Fig. 8. There is a strong similarity between this test case and 
the one shown in Fig. 6. This resemblance and also comparing the form of Ct,±i/ 2 , Eq. (4.3), with the new FV 
smoothness indicator T ], Eq. (3.26), strongly suggests to take a closer look at £. In recent years, this parameter has 
attracted a lot of attention, see e.g. [30], [1], [2], [15] and references therein. Originally, £ was introduced by Jiang 
and Shu to avoid the denominator to become zero [12]. The authors called it a small positive number and set to 
£ = 10 6 for their test-case studies. Therefore, WEN03 with the fixed value £ = 10 6 will be called WENO-JS 
henceforth. There are several drawbacks by fixing £ to some value Eq. One of them is the following scenario which 
might occur since 8 , i = &(Ax) in smooth parts of the solution: 

1. For large grid sizes, 5.+ 1 | 2> £o holds, leading to the homogeneous function 7 /weno3 —► ,, 5 i), 

Eq. (4.6), and thus to low order. 

2. Refining the grid leads to 1 8 j± 1 1 <C £o, and yields W\vhno 3 —> H. (8 i i, 8 j+ 1 ), Eq. (4.5), which is the full-3 rd 
order reconstruction. 

This phenomenon is also demonstrated in Fig. 8. Wanting £ to be adaptive for all grid discretizations means that it 
has to depend on the grid size, i.e. £ = e(Ax). This is also what Yamaleev and Carpenter claim in [29,30]. However, 
the authors do not simply replace £ by a function depending on Ax, they rather define new weight functions 

(Xi ( X \ 

w k =—I —, <4 = Yk\ 1 + — 7T , fc = 0,...,r-l, (4.7a) 

I- = o«, V £ + PkJ 

£ = max(11111? II( m o) 2 ||ij ■ ■ •, ll(«o^ 1) ) 2 ||i) Ax2 > ( 4 - 7b ) 

where £2 f / is a set of points where the initial condition is discontinuous. In this paper we consider r = 2. Thus, 
yo = Yi - i / 2 - 7 i = Yi+ 1 / 2 > etc. are defined as before and x is the square of the undivided difference on the entire 
stencil. In case of the compact three-point stencil this is x = (h,+i — 2i7, ■ + m,-i ) 2 = — A - 1 ) 2 - N° te tb at 
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the definition of e is not invariant under translation of the initial condition uq. This leads to a different limiting 
behavior even though no may simply be shifted by a constant, cf. Sec. 6.3. 

Arandiga et al. [2] show that using the weights proposed by Yamaleev and Carpenter, the resulting scheme reaches 
the maximal order (2 r— 1) for sufficiently smooth solutions. However, near discontinuities, the scheme achieves 
order of accuracy Ax 2 ) which is worse or equal to 8 )’(Ax r ). r > 2, the order of accuracy of the underlying ENO 
scheme. Arandiga et al. report to have fixed this issue by slightly changing the weight functions using 


“ , = b ( 1 + Gt&) )■ k=0 ." = 

£ = KAx q , with K > 0, q € N, q < 4r — 4 — rfol. 


(4.8a) 

(4.8b) 


In the case of the three-point stencil, fx = 1 and therefore the weight formulation remains the same. In the numerical 
test cases carried out in [2], K is set to 1, which makes e a dimensional quantity that might be affected by changes 
of units. We will not consider these weight functions in our numerical experiments in Sec. 6 since the definition 
of e in Eq. (4.7b) is clearly more elaborate. Solely in Sec. 6.4 we compare the obtained results with the scheme 
setting K— 1, i.e. e = Ax 2 to show the importance of the definition of e. 


5 A Unifying View 

In this section we want to place the different methods in a unifying view to point out their differences and especially 
their similarities. 

For the sake of simplicity, we do not discuss the three-dimensional plots of the limiter functions but rather sectional 
views at fixed values of 5 ; i. This means, the limiter functions are depicted with a one-dimensional dependence 
on 5 ( _ i. Lemma 3.2 states that for 5 1+ 1 = 1 this is equivalent to the standard (0, 0(0))-plots such as Fig. 2, and 
found in textbooks. 

To give an overview. Fig. 9 shows the limiter functions treated in the previous sections for 5. + i £ {2,1,0.5,0.1}. 
In Fig. 9 we can nicely see that all limiter functions satisfy Property 3.1 (i). This is, for 0 = 1, i.e. 8 . i = 5 i, the 
limiter functions continuously go through the point (5 ;+ i .5 i), which corresponds to 0(0) = 1. This situation 
represents smooth parts of the solution. It can be seen that WENO-JS has slope j only close to the point 5. | = 
8 : , i meaning that as soon as the slopes (of same sign) differ from each other, the limiting takes effect, even though 
the function is still smooth. The smaller the values of (5 ( - ±1 / 2 , the stronger is this effect, see the trend from Fig. 9a 
to 9d. 

For 8 i —r ±°°, i.e. 0 —> ±°°, all considered functions, except H \, tend to a constant. For the WENO schemes, 
this constant is non-zero, contrary to the negative part of the FV limiter. However, for the reconstruction of the 
cell interface values, H{8 j i, 5. i) = const, leads to first-order accuracy, independent of the value of the constant. 

In Fig. 9c and 9d we can see the switch of the FV limiter to the full-3 rd order reconstruction, highlighted by the 
dashed black line. Clearly, the construction of the transition zone is rather ad-hoc for the FV limiter, while it comes 
more naturally in WEN03. However, the FV setup allows a much more systematic control about the shape of the 
function H far away from the origin, particularly in the regions where <5. i and 5 i have opposite signs. It is 
here that the FV and WENO limiters show a very different behavior. In the point 8 j _i = 8 /+ i, i.e. at extrema, 

all proposed limiters equal Hj,. However, while //(}■’ has the same slope as the full-3 rd order reconstruction, the 
WENO limiters have negative slopes. This indicates that the FV limiter reconstructs extrema with higher accuracy, 
since in general extremal points might not lie exactly at 8 . i = — 5., i but rather 8 i « 8 i. In these cases, 

the new limiter function approximates the full-3 rd -order reconstruction, leading to higher order solutions. 

Another interesting observation is the nature of and WENO-YC at 5._i = 0. When considering the univariate 
limiter functions 0(0) with 0 = 5 ; _ i /8 j+ 1 , the point 5 i = 0 is equivalent to 0 = 0. In this case, the conventional 
FV limiter functions in the MUSCL framework are set to 0, yielding a first-order reconstruction. This is also the 
case for WENO-JS, whereas the limiters proposed in this work and by Yamaleev and Carpenter yield non-zero 
contributions for small 8 . i, see Fig. 9b, 9c, and 9d. This is exactly the desired behavior close to extrema where 
one slope may be zero and the other slope small but non-zero. Further away from the origin in the (5 ; _i ,5 (+ i)- 
plane, it is clear that the limiter functions should yield 0 for 5._ i = 0 because this situation may correspond to a 
discontinuity. This feature can be observed in Fig. 9a. 
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(c) 5,,i = 0.5. (d)5, i=0.1. 

1 ' 2 2 

Figure 9 - Sectional view of different limiter functions for the fixed values 5 ;+ i € {2,1,0.5,0.1}. 


6 Numerical Results 

In this section, we apply and compare the different limiter functions discussed in the previous sections. In all test 
problems proposed in this work, we compare 

1. the original WEN03 scheme as introduced in [12], i.e. fixing e = 10 6 , called WENO-JS; 

2. WEN03 with the weight functions proposed by Yamaleev and Carpenter [30], called WENO-YC; 

3. the full 3 lcl -order reconstruction # 3 , Eq. (3.11); and 

4. the FV limiter function //} L , as introduced in Sec. 3. 

The time derivative is approximated using the 3 ld -order Runge-Kutta method as described in [24]. 

For all test cases, rather than presenting tables with errors and convergence rates, we plot the L\ - and L„-errors. 
In these plots, we include reference slopes of the design order of accuracy, i.e. the order of accuracy the schemes 
obtain in theory. 

6.1 Preliminary Test Case 

(c) 

In this test, we want to point out the importance of the choice of e for WENO-YC as we have done for H C j and 
WENO-JS in Sec. 3.2 and 4.2, respectively. We set e =CAx 2 . According to Eq. (4.7), C = max^o rf (||nQ||, ||(«()) 2 ||). 
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-*■ e=1000.0 Ax 2 
-e- e=1.0 Ax 2 
e=0.1 Ax 2 
-g- 6=0.01 Ax 2 
6=0.001 Ax 2 


Figure 10 - Zoom of the solution of the preliminary test case, Eq. (6.1). This has been obtained with the 
WENO-YC scheme, £ = CAx 2 , C G {10 3 ,1.0,10 10 2 ,10 3 } with n = 80 grid cells and 

^end = 1-0, 


However, in this test case we set C G {10 3 ,1.0,10 ',10 2 ,10 1 } to study the influence of the coefficient. The re¬ 
sults for the repetition of the test 


u t + u x = 0 (6.1a) 

uo(x) = sin(^jr), x€ [-1,1] (6.1b) 

with n = 80 cells and f e nd = 1.0 is depicted in Fig. 10 for different values of £. One can clearly see the loss in 
accuracy if the constant for £ is chosen too small. Note that the correct value for £, as proposed by Yamaleev 
and Carpenter [30], is £ = max(l, n 2 ) Ax 2 = n 2 Ax 2 . This indicates that for smooth functions a coefficient which 
is too small decreases the quality of the approximation. This can also be observed in Fig. 11 which shows the 
L\- and L,„-errors of the solution. We can see that the solution with £ = 1000 Ax 2 is 3 rd -order-accurate in both 
norms starting at the lowest resolution. While the solution with £=1.0 Ax 2 is still directly 3 rd -order-accurate in 
the L |-norm starting from n = 20 grid cells, we observe that the smaller the coefficient of £, the larger n must be 
chosen so that 3 rd order is achieved. 

(c) 

Of course, the same conclusion holds true for the choice of a in This test thus shows that we need to carefully 
evaluate Eq. (4.7) and (3.26). A misinterpretation of the coefficients may lead to results which are significantly 
worse than what the scheme is capable to achieve. 




n=tt cells n=tt cells 

Figure 11 - Double logarithmic plot of the L\- and L„ -errors versus number of grid cells of the solution 
obtained with 7/wenoyc- advected until f en( j = 1 with CFL number v = 0.9. Numerical 
solution for different values of £ = CAx 2 , C G {10 3 , 1.0, 10 1 .10 2 .10 3 }. Left: Li-error, 
right: -error. 
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(a) Solution of Eq. (6.2), (6.3) with different 
numerical schemes, slightly zoomed in. 



(b) Enlarged view of the maximum of the solution. 


Figure 12 - Results of advection equation Eq. (6.2) with smooth initial condition, v = 0.8, n = 170 grid 
cells and f en( j = 10, i.e. the solution has been advected 10 times around the domain. 


6.2 Advection Equation with Smooth Initial Condition 

With this first test case we aim to verify that all considered schemes are 3 ld -order accurate for smooth solutions. 
We solve the linear advection equation, 

u t +u x =0 ( 6 . 2 ) 


with the smooth initial condition 

, f (0.5 + 0.5cos(5;r(x-0.5))) 4 if 0.3<x<0.7 

«w = { 0 else - - < 63 > 

and periodic boundary conditions. The computational domain is [0,1], the CFL number v = 0.8 and the so¬ 
lution is advected until f en( j = 10. The spatial resolution is the sequence of refined uniform grids with n = 
20,40,50,100,120,170,200,300,500,700,1000,1500,3000 cells. For WENO-YC, according to Eq. (4.7), we set 
£ = 20.67 Ax 2 with Ax = 1 /n. For the FV limiter function Eq. (3.28) and (3.26) we fix a = max A . e£1 \Q rf |mq(x)| = 
493.48. 

Fig. 12a and 12b show zooms of the solution at f e nd with different numerical schemes on a 170-cell grid. It can 
be seen that and WENO-YC, as well as the full-3 -order scheme, perform much better than the conventional 
WENO-JS scheme in terms of accuracy. This can also be observed in Fig. 13 which shows the L\- and -errors 
obtained at f e nd- The convergence rates of our proposed FV scheme, WENO-YC, and the unlimited 3 ld -order 



Figure 13 - Double-logarithmic plot of the L\ - and Loo -errors vs. number of grid cells for the advection 
equation with smooth initial condition, Eq. (6.2), (6.3). 
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Figure 14 - Enlarged view of the discontinuities of the solution of Eq. (6.2), (6.4) with different numerical 
schemes using n = 320 cells, v = 0.8,f e nd = 10. 


scheme reach 3 rd order starting at n = 150 grid cells whereas WENO-JS shows this order of convergence only for 
very large numbers of cells and with a larger error constant. 


6.3 Advection Equation with Discontinuous Initial Condition 

In this section we want to discuss the behavior of the numerical schemes for solutions containing discontinuities. 
Therefore, we consider the advection equation, Eq. (6.2) with a = 1 and square wave initial condition 


u 0 (x) 


1 for — 0.5 <x < 0.5 
0 else. 


(6.4) 


The computational domain is [—1,1], the CFL number v = 0.8 and the solution is advected until r en d = 10, which 
corresponds to 5 periods in time. Due to the large gradients contained in the initial condition, the limiter functions 
have to take effect in order to avoid spurious oscillations to appear. Solving this test case with the full-3 ld -order 
reconstruction Eq. (3.6) generates oscillations. This is the reason the FV limiter functions presented in Sec. 3 
restrict the reconstruction to l st -order in these situations and WEN03 reduces to Eq. (4.6). The WENO-YC 
parameter e is given by Ax 2 (see Eq. (4.7b)) and reduces to Hu because a = max^m^ Mq (x) = 0. 

This test case nicely shows the already mentioned drawback of the definition of e in WENO-YC, Eq. (4.7b), 
namely the coefficient of e which is not translationally invariant if the initial condition is shifted. To point this out, 
a second test case has been chosen, where the initial condition has simply been shifted by +100, i.e 


IC+ioo : Mo.+iooW 


101 for —0.5 < x < 0.5 
100 else. 


(6.5) 


When applying WENO-JS to this new initial condition, the solution is does not show more oscillations, because 
£ is fixed to 10 6 . In fact, the solution is the same as for liq(x), only shifted by +100. Also, in the proposed FV 
limiter, the constant a = max^m^ I m o+ioo( x )I = ® does not c h an 8 e - Thus, the scheme yields the exact same 
results, shifted by +100. However, for WENO-YC, £, as given by Eq. (4.7b), is no longer Ax 2 but 20201 Ax 2 . 
The higher value of £ leads to augmented oscillations in the solution, as can be seen in Fig. 14. Here, for better 
comparison, the solution w+iooCx+end) °f the test with shifted initial condition, which lies in the range [99.9,101.1], 
has been shifted to the range [—0.1,1.1]. Thus, both test cases, Eq. (6.4) and (6.5) lie between —0.1 and 1.1 in the 
plots and the magnitude of the oscillations can be compared. As seen in Fig. 14, the original WEN03 scheme does 
not cause any oscillation but it is rather dissipative. Our proposed FV limiter does not produce overshoots either. 
Additionally, it approximates the sharp structure of the initial condition better than WENO-JS. WENO-YC with 
correctly chosen parameter approximates the discontinuity almost as good as J/ 3L , however, it causes oscillations. 
These are even larger for the badly chosen £. This behavior can also be observed in Fig. 15a and 15b which show 
the L \-error and the total variation (TV), respectively. The best order of convergence we can expect from a 3 rd - 
order scheme and a solution containing discontinuities is 3 /4. This can be shown by the Fourier analysis of the 
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(a) Double-logarithmic plot of the L \-error vs. number 
of grid cells. 


(b) Total Variation vs. number of grid cells of the 
different schemes. 


Figure 15 - Results of advection equation Eq. (6.2) with discontinuous initial condition., v = 0.8 and 
tend =10, i.e. the solution has been advected 10 around the domain. 


modified equation with self-similar initial condition, cf. [17], 

Even though, both solutions computed with WENO-YC cause oscillations, they obtain the order of accuracy 3/4. 
The more dissipative WENO-JS scheme is also of order 3/4 but with a larger error constant. Among the tested 
schemes, the best error constant is obtained with our proposed FV limiter function. The total variation of all 
schemes represents their behavior as seen in the solution. WENO-JS attains the TV of the exact solution TV (u ex ) 
from below, meaning that is does not cause overshoots at all, whereas WENO-YC is larger than TV(u ex ) for all 
spatial discretizations. is closest to TV(u ex ) and lies never above TV(u ex ) for all time steps, i.e. it does not 
cause oscillations at any time. 

6.4 Initial Condition Containing Smooth and Non-Smooth Features 

We consider the same setup as in Sec. 6.2 with CFL number v = 0.8, f e nd = 10.0 and initial condition 

M 0 (x) =max(min(^/ [ - -2,- -2^) +l,o) +exp j sin(307Tx). (6.6) 

In this problem, we want to test how accurate the different schemes resolve small features in a larger setting of a 
more complex solution. The spatial discretizations are n = 20 • 2‘, i = 0.... ,7 grid cells. For WENO-YC, we set 
£ = 1042.83 Ax 2 , Ax = 1 /n, according to Eq. (4.7b). For the FV limiter function we fix a = ma |wq(x)| = 
8887.87. We run a third case with WENO-YC but setting e = Ax 2 to show the difference in performance as 
described in Sec. 4.2 and 6.1. This test case corresponds to the weight functions proposed in [2, 15], where 



Figure 16 - Initial condition (6.6) containing smooth and non-smooth features. 
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-a- WENO-JS 
-<© WENO-YC, 
e=1042.83 Ax 2 
■+- WENO-YC, 
e=Ax 2 

-e-full 3 rd order 

Uex 

X X 

Figure 17 - Results of different schemes for the advection equation with more elaborate initial condition: 

Zoom of two significant regions of the solution with n = 640 grid cells at f en d = 10, v = 0.8. 




e = KAx 2 was used with K = 1. Furthermore, a test with the full-3 rd -order reconstruction f/ 3 , Eq. (3.11), was 
performed and compared to the numerical schemes discussed in this paper. 

Fig. 17 shows a close up view of two significant regions of the solution with n = 640 grid cells. We observe 
that our proposed FV limiter function and WENO-YC with correctly-chosen parameter e perform much better 
than WENO-YC with e = Ax 2 . The results of these two schemes are very close to the ones of the full-3 rd order 
reconstruction. Also, all four schemes outperform the conventional WENO-JS scheme. This can also be seen in 
Fig. 19, which shows the double-logarithmic plot of the L\ - and L„,-errors versus the number of grid cells in the 
smooth part of the solution, x £ [0.4,1], The solution cannot be accurately represented with few grid cells by any 
of the treated schemes. Even the full 3 rd -order reconstruction does not resolve the details and therefore has a large 
error constant. As soon as a reasonable space discretization is reached, the order of convergence reaches 3 rd order 
if only the range x £ [0.4,1], i.e. the smooth part of the solution, is regarded. Solely WENO-JS does not reach this 
order even at the highest resolution. 

If the error is considered on the whole domain [0,1], the convergence rate of all schemes degenerates for higher 
resolutions. This is due to the fact that at higher resolutions the kinks in x £ [0.1,0.3] become well-resolved and 
thus recognized as non-smooth. The initial conditions are an interesting test case combining smooth and non¬ 
smooth features and therefore testing the capabilities of limiter functions. As shown in Fig. 18, near the position 
of the kinks, all schemes - with the exception of WENO-JS and H 3 l - generate undershoots. For our proposed FV 

(c) 

limiter this can be explained with the large asymptotic region. Since the maximal second derivative of the 
initial condition is very large, a = max A . en \Q rf Mq (x) = 8887.87, the region where reconstructs with full-3 rd 
order is relatively large. At the same time, the discrete second derivative of the kinks of the triangle are small 
compared to the extreme regions of the smooth part, and at these points, the solution is reconstructed using // 3 
with no limiting. As a result, //) L causes undershoots, just as // 3 does, because this is what is efectively used at 
the kinks. However, compared to other limiting methods, we can in principle control these undershoots using our 



~H 3L 

& 

WENO-JS 
-9- WENO-YC, 
6=1042.83 Ax 2 
WENO-YC, 
e=Ax 2 

-©-full 3 rd order 

Ucx 


Figure 18 - Zoom of the kink of the solution with n = 640 grid cells at f e nd = 10, V = 0.8. 
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Figure 19 - Double-logarithmic plot of the L\ - and Lc-norm of the error vs. number of grid cells in the 
range x £ [0.4,1], i.e. the smooth part of the solution of (6.2), (6.6). 


(c) 

FV limiter . By choosing a smaller value for a, i.e. a smaller asymptotic region, these undershoots can be 
completely avoided, as can be seen in the test case with pure discontinuity, cf. Sec. 6.3. An adaptive value for a 
would therefore eliminate the undershoots in the non-smooth region x £ [0,0.4] while still resolving the smooth 
smooth parts in x £ [0.4,1] with high order accuracy. However, such a local adaptivity would necessarily require 
to use a wider stencil for the reconstruction, because as shown in Fig. 5, three points can not distinguish between 
the kink and a strongly curved extremum on a coarse grid. 

6.5 Sod Shock Tube Problem 

Let us consider Sod’s problem, which describes a shock tube containing two different ideal gases at the left and 
right side of a diaphragm, placed at x = 0. The density, velocity, and pressure of the gases in the left and right 
region are given by 

Pl\ ( l-0\ (p R \ / 0.125,\ 

vl = 0.0 , v* = 0.0 . (6.7) 

Pl) \1 -0/ \Pr) V 0.1 / 

At time t > 0, the diaphragm is removed and the gases begin to mix. The time evolution is described by the one 
dimensional Euler equations, 

u,+f(u) v = 0 (6.8a) 

with u=(p,pu,E), the flux function 

f(u) = (pu,pu 2 +p, u(E+ p)) T , (6.8b) 


and the equation of state for ideal gases 

E=-l—- + -pu~, (6.8c) 

y— 1 2 

with y = 1.4. The computational domain is set to [—2,2] and the test is conducted with N = 100 grid cells 
until Tend = 0.8. We compare the new limiter function //,, with the full 3 -order reconstruction and WENO-JS. 
WENO-YC is not included in the plots because it produced negative values for pressure, when run with e = 2.25. 
Due to the purely discontinuous form of the initial condition, we obtain a = 0, just as in Sec. 6.3. This means, that 

(c) 

Z/j L = fp[, i.e. no asymptotic region exists in this test case. The reconstruction techniques have all been applied 
to the primitive variables p, u and p. 

Sod’s shock tube problem leads to three characteristic waves, which can be seen in the solution. Fig. 20. Both, the 
density and pressure profile show the rarefaction wave and the shock. The contact discontinuity can only be seen in 
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X X 

(a) Density profile. (b) Pressure profile. 


Figure 20 - Solution of different reconstruction techniques for Sod’s Problem, Eq. (6.7), (6.8) on the 
domain [—2,2] with N = 100 grid cells, CFL v = 0.95 at r enc | = 0.8. 


the density profile. Fig. 20a, not in the pressure distribution. Fig. 20b. The solution shows that applying the full-3 rd 
order reconstruction leads to over- and undershoots close to the discontinuities. WENO-JS yields good results 

(c) 

concerning this feature, however, does not approximate the true solution as accurate as L , see also discussion in 
Sec. 6.3. 


6.6 Shu-Osher Problem 

In this problem, originally introduced by Shu and Osher [25], a Mach 3 shock is interacting with sine waves in the 
density profile. The computational domain is fixed to [—4.5,4.5] and the shock at time t = 0 is situated at x = —4. 
The initial conditions of the primitive variables density, velocity and pressure, to the left and right of x = —4 are 
given by 



x 


Figure 21 - Solution of different reconstruction techniques for the Shu-Osher Problem, Eq. (6.9), (6.8) with 
N = 640 grid cells on the domain [—4.5,4.5], CFF v = 0.95 at time Tend = 1.8. 
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Figure 22 - Solution of different reconstruction techniques for the Shu-Osher Problem, Eq. (6.9), (6.8)with 
N = 1280 grid cells on the domain [—4.5,4.5], CFL v = 0.95 at time 7^ = 1.8. 


(Pl\ /3.857143\ ( p R \ (\ +0.2sin(5x),\ 

\v L = 2.629369 , v* = 0.0 , (6.9) 

\PlJ \10 33333 J \p R J V 10 / 

and the time evolution is governed by the one-dimensional Euler equations, Eq. (6.8). Just as in [25], the solution 
is computed at 7 enc j = 1.8. We compare the solutions of the full 3 rd -order reconstruction, WENO-JS, WENO-YC 

(c) 

(where e = 21.932), and the new limiter function // , L with N = 640 and N = 1280 cells to a reference solution, 
which is the numerical solution of WENO-JS with 10,000 grid cells. The CFL number is set to v = 0.95 in all 
tests. 

The combined limiter function includes an asymptotic region, cf. (3.26), with a = 5.0 the reconstruction 
techniques have all been applied to the primitive variables p,u and p. 

Fig. 21 shows a comparison of the solutions obtained on 640 grid cells of the full 3 ld -order reconstruction, WENO- 

(c) 

JS, WENO-YC, and the new limiter function 7/ 3L . A zoom of the areas of interest of the solutions computed on 
1280 grid cells is shown in Fig. 22. Overall, it can be seen that applying the full-3 rd order reconstruction or WENO- 
YC leads to under- and overshoots close to discontinuities. WENO-JS does not produce overshoots, however, it 
limits too much, so that the reference solution is not approximated as well as by the limiter function 7/ 3L . This is 
especially visible in regions with large gradients. 

7 Conclusions 

In this paper we have analyzed 3 rd -order finite volume and WENO schemes. These schemes reconstruct the so¬ 
lution at cell interfaces, each reconstruction based on only three mean values. We have then placed the different 
schemes in a unifying context. More specifically, we have analyzed and improved the FV limiter by Cada and 
Torrilhon [4]. Our proposed FV limiter does not contain artificial parameters anymore. For a unifying view, 
we considered the 3 rd -order WENO schemes proposed by Jiang and Shu [12] and Yamaleev and Carpenter [30]. 
Plotting all schemes in the (5._i, 8. i i-plane revealed certain similarities which could also be observed in the 
numerical experiments. For smooth solutions, the FV limiter and WENO-YC show equally good results. How¬ 
ever, near discontinuities, the numerical experiments showed some oscillations using the WENO-YC scheme, as 
predicted by Arandiga et al. [2]. This effect is avoided using our proposed FV limiter function, which reduces 
oscillations near discontinuities significantly. 
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